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ABSTRACT 

A linear algebraic solution is provided for the problem of retrieving the location and time of occurrence of 
lightning ground strikes from an Advanced Lightning Direction Finder (ALDF) network. The ALDF network 
measures field strength, magnetic bearing, and arrival time of lightning radio emissions. Solutions for the plane 
(i.e., no earth curvature) are provided that implement all of these measurements. The accuracy of the retrieval 
method is tested using computer-simulated datasets, and the relative influence of bearing and arrival time data 
on the outcome of the final solution is formally demonstrated. The algorithm is sufficiently accurate to validate 
NASA's Optical Transient Detector and Lightning Imaging Sensor. A quadratic planar solution that is useful 
when only three arrival time measurements are available is also introduced. The algebra of the quadratic root 
results arc examined in detail to clarify what portions of the analysis region lead to fundamental ambiguities in 
source location. Complex root results are shown to be associated with the presence of measurement errors when 
the lightning source lies near an outer sensor baseline of the ALDF network. For arbitrary noncollinear netw r ork 
geometries and in the absence of measurement errors, it is shown that the two quadratic roots are equivalent 
(no source location ambiguity) on the outer sensor baselines. The accuracy of the quadratic planar method is 
tested w r ith computer-generated datasets, and the results are generally better than those obtained from the three- 
station linear planar method when bearing errors are about 2°. 


1. Introduction 

Advanced Lightning Direction Finder (ALDF) sen- 
sors, developed by Global Atmospherics Inc. (GAI), 
have the ability to detect the field strength, magnetic 
bearing, and arrival time of lightning radio emissions. 
In 1992, Lightning Location and Protection, Inc. (a di- 
vision of GAI) completed development of an Improved 
Performance from Combined Technology (IMPACT) 
method for determining the location and time of oc- 
currence of lightning return strokes from these data 
(Cummins et al. 1993). The IMPACT algorithm is based 
on minimizing a x 2 function similar to that provided in 
Eq. (1) of Hiscox et al. (1984) but generalized to ac- 
commodate arrival time data. The lightning time of oc- 
currence /, the longitude A, and the latitude <p of the 
lightning source on an ellipsoidal earth is estimated. 

Since the IMPACT algorithm uses a numerical ap- 
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proach to determine the absolute minimum of the non- 
linear x 2 hypersurface, it does not represent an analytic 
solution to the problem; that is, the source location and 
time of occurrence are not directly determined in terms 
of the measurements and measuring network geometry. 
Instead, the algorithm uses the computer to search for 
the optimum values of (A, <p, t) that minimize x 2 \ each 
new set of measurements implies starting an entirely 
new search for an answer. In effect, the actual solution 
is estimated using the constraints of the data and the 
power of the computer. 

Relative multiple minima in the ^ 2 (A, <p, t) hypersur- 
face can lead to a premature termination in the computer 
search and an erroneous solution. The presence of data 
errors can generate additional relative minima and ad- 
ditional solution errors. Since it is computationally ex- 
pensive to check for erroneous solutions during real- 
time processing of ALDF data, lightning source solu- 
tions depend, in general, on where in the solution space 
(A, <p, t) the computer search begins. 

Although the IMPACT algorithm has had successful 
practical application (Cummings et al. 1998), the spe- 
cific algorithm software is proprietary and is not widely 
distributed free of charge to the scientific community. 
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In this study, our interest was to develop an economical, 
yet useful (four-station) ALDF ground-truth site at the 
Melville Island-Darwin, Australia, region — and else- 
where in the world — that could be used to validate the 
National Aeronautics and Space Administration’s 
(NASA’s) space-based lighting detectors [the Lightning 
Imaging Sensor (LIS), and the Optical Transient De- 
tector (OTD) described in Christian et al. (1992) and 
Goodman et al. (1995), respectively!. We desired a com- 
putationally quick algorithm that easily ingests all types 
of ALDF measurements, and that produces ground- 
strike locations with accuracies better than 4 km (the 
nadir resolution of LIS) for sources that are within a 
few hundred kilometers of the ALDF network. This has 
been the primary motivation behind our algorithm de- 
velopment. As such, we did not require an algorithm 
that accounts for earth curvature. 

A secondary motivation behind this work was real- 
ized during our algorithm development. We learned that 
a serious pedagogical effort to examine the mathemat- 
ical foundations of lightning ground-strike retrievals 
does not appear in the literature. Perhaps the wide use 
and practicality of powerful numerical computer algo- 
rithms have offset the desire to obtain analytic solutions. 
In retrospect, the proper mathematical approach is to 
first solve the least difficult form of the retrieval problem 
in a thorough way, and then methodically upgrade the 
formalism to account for additional physical complex- 
ities such as associated with the addition of new data- 
sets, earth curvature, and wave propagation effects. An 
effort must also be made to clearly explain fundamental 
ambiguities that arise in the retrieval process when an 
experimenter has limited measurements available due 
to sensor failure or other causes. Although there have 
been some analytic efforts to account for earth sphericity 
for magnetic bearing data inversions (e.g., see Orville 
1987) a complete pedagogy that ultimately leads to rig- 
orous retrievals of ground strikes on an ellipsoidal earth 
surface (and that uses all ALDF data constraints in a 
coherent/simultaneous and optimum fashion) is notably 
absent in the literature. 

To the best of the authors’ knowledge, no one has 
yet explicitly provided an analytic solution to the prob- 
lem of determining lightning source location and time 
of occurrence using collective measurements of field, 
magnetic bearing, and arrival time measurements when 
earth curvature and propagation effects are neglected. 
This is a fundamental starting point. For multiple dataset 
inversions, the relative importance of bearing and arrival 
time data on the outcome of the final solution has not 
been formally demonstrated. Furthermore, no thorough 
investigation of solution ambiguities have been provid- 
ed when one is limited to just three arrival time mea- 
surements. This writing introduces theoretical deriva- 
tions that address each of these problems. In so doing, 
we arrive at the required OTD-LIS ground truth algo- 
rithms. 

We specifically determine the source location (x, y) 


as a mathematical function of the measurements under 
a variety of conditions (i.e., differences in the number, 
location, and type of measurements). Since algebraic 
solutions are obtained, we do not need to invoke a com- 
puter search algorithm to determine optimum solution 
parameters. We provide unique physical insight into the 
nature of the retrieval problem because we determine * 
exactly how the measurements are specifically related 
to the lightning source location (and time of occurrence). 

A Linear Planar (LP) method is first introduced that 
allows one to simultaneously analyze field, bearing, and 
arrival time measurements, and a means for optimally 
weighing the bearing data relative to arrival time data 
is demonstrated. To the best of our knowledge, the linear 
approach provided here has not been attempted else- 
where in the literature. The method involves one large 
system of linear equations that offers a high degree of 
flexibility from the point of view of the user’s appli- 
cational needs. For example, if only a certain number 
and type of measurements are available in an experi- 
ment, the linear system of equations degenerates into a 
smaller set of equations, and a straightforward solution 
process is retained. 

We also introduce a Quadratic Planar (QP) method 
that can be used when only three arrival time measure- 
ments are available. Such a situation arises if there are 
sensor hardware failures and/or when field amplitude 
and bearing measurement data quality is unacceptable. 
Although this method is mathematically nonlinear, full 
analytic solutions are derived. Physical insight about the 
nonlinear solution space, not discernible from conven- 
tional x 2 analyses, is fully described by examining in 
detail all quadratic root solutions derived from the QP 
method. For example, we show explicitly that a certain 
mathematical discriminant vanishes for certain lightning 
locations, and that these source locations produce com- 
plex roots (negative discriminants) in the presence of 
measurement errors. 

Extensive tests of the LP and QP retrieval methods 
are provided using computer-simulated datasets and 
these methods are applied in a study of ALDF data that 
were obtained from the Maritime Continent Thunder- 
storm Experiment (MCTEX) analysis region in Darwin, 
Australia (Keenan et al. 1994, 1996). Data from this 
network comprise one of several ground-truth sites for 
the validation of OTD and LIS. 


We begin by considering n > 3 sensors situated at 
locations r,, / = 1, 2, . , . , n relative to some origin. 
Each sensor has the capability to measure the arrival 
time f., magnetic bearing and field strength F, of 
the radio emissions from a lightning source with loca- 
tion, r, time of occurrence, /, and radiation source 
strength, s. Hence, from the 3 n measurements {(f 1# <f > ,, 
F,X . . . , (C, 4>„, F„)} we wish to determine the five 


2. Linear Planar (LP) method 
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Pig. 1. Geometry associated with the LP method. 


unknowns (x, y, z y t, s). In so doing, we neglect earth 
curvature. 

Figure 1 summarizes the geometry of the LP model. 
Because ALDF sensors might not be deployed on a flat 
topography, the i th sensor located at r ( need not lie in 
the xy plane, that is, z- t ^ 0 in general. The relative 
position vector follows standard physics convention, 
that is, it points from the source at r to the observation 
point r,, so that R, = r, - r. Neglecting refractive effects 
in the atmosphere, the excitation time of the /th sensor 
is 

t, = t + -tf,., (1) 

c 

where c is the speed of light. Solving for tf,, squaring, 
and rearranging terms leads to 

“(rf - c 2 t j) = x t x + y\y + z 4 z ~ c 2 t,t 

- - c 2 ry ( 2 ) 

It is desirable to remove the last term on the right-hand 


side of (2) since it is nonlinear in the space and time 
variables. To do this, we define the measurement 

«, - - c'-t}) - |(rr - c'-t r). (3) 

A comparison of (2) and (3) shows that a , is linearly 
related to the lightning location, r = (x, y, z ), and light- 
ning time of occurrence, t\ that is, 

a 4 - (x, - xjr + (y, - }>,)>' + (z, ~ z x )z 

- c 2 (t, - t,)t\ i = 2, 3, ... , il (4) 

A detailed investigation of this linear form has been 
provided in Koshak and Solakiewicz (1996). 

Next, we consider the information content of ALDF 
bearing data. From Fig. 1 we see that the lightning lo- 
cation (x, y) is given by 

x = Xj + p t cos (f), y = y t + p t sin</>,, (5) 

where p, is the horizontal distance from the /th site to 
the lightning ground-strike location. It is useful to define 
the measurement 

P, = x, sin</>, - y t cos <p r (6) 

Solving (5) for x, and y ( and substituting into (6) gives 

P, = (sin <j> t )x - (cos <f>,)y. (7) 

Finally, we consider measurements of the radiated 
field strength. Assuming a 1/tf, attenuation in the ra- 
diation field gives 

F, = ^ (8) 

Once again, we solve for tf,, square, rearrange terms, 
and define the measurement 

y, = \ir} - rf). (9) 

This leads to the following relation: 

y, = (x, ~ ■*,)* + (>■, - >’,)>' + (z, - z,)z 

* i{j} ' ^ ' = <l0) 

If we consider only n — 3 sensors, (4), (7), and (10) 
can be combined to give 
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where 8 is a weighting factor chosen as 10" m, 'Jc = 
5 [(F,/F,) 2 — 1] is a dimensionless parameter, d t = ct, 
d, = and £ = [F*(x 2 — X|) 5/2 ] 1 is a scaling factor 
to be described below. Defining the column vector on 
the left-hand side of (1 1) as g, the matrix by K, and the 
remaining vector by f, we may rewrite (11) as 

g = Kf. (12) 

All elements of K and f are in units of meters, and 
all elements of g are in squared meters. This was ac- 
complished by retaining a factor of c in front of the 
time-difference measurements in the first two rows 
(fourth column) of K, by multiplying (7) by the weight- 
ing factor 8, and by scaling the field F, and source 
strength s, each by the factor that is, by making the 
substitutions F, — » £F,, and s — » £s in (8). (Note that 
the sites must be numbered in such a way that * 2 > x, 
so £ is not complex; this can always be accomplished 
since the numbering of sites is arbitrary and because 
the translation and rotation of the x-y coordinate system 
used in the LP method is arbitrary.) 

In general, the K matrix has 3n — 2 rows and 5 
columns, where n — 1, 2, 3, .... If there is only n = 
1 sensors in the network, K degenerates into a row vector 
and (12) is underdetermined. If there are n = 2 sensors, 
K will have 4 rows and (12) will still be underdeter- 
mined. For n > 3 sites, K will have >7 rows and (12) 
will be overdetermined. For overdetermined systems, f 
can be retrieved using the least squares inversion pro- 
vided in Twomey (1977): 

f = (KK) ‘Kg, (13) 

where the tilde represents matrix transposition. The 
source time of occurrence and source strength are de- 
termined as t ~ djc, s = d\ f2 fi;, respectively. 

From the foregoing generalities, we now note that 
ALDF sensors trigger on the initial upward current surge 
that occurs close to the ground level in the return stroke 
[see, e.g., the transmission line model discussed in 
Uman et al. (1973)]. Hence, the source can be regarded 
as being located at z = 0. In this case, we can remove 
the third component of f, that is, we consider the column 
vector f = col(x, y, d n d s ), and we remove the third 
column of K. We then regard the expression in (12) as 
a (3n — 2) by four system of linear equations. In this 
case, n = 2 sensors generates a (4 X 4) K matrix so 
that (12) is a determined system with direct solution f 
= K g. Hence, the LP method can be used by an ex- 
perimenter that has only two sensors, each measuring 
bearing, arrival time, and field amplitude. In this case, 
source location (x, v), time of occurrence (, and source 
strength s can be retrieved. If the two sensors do not 
provide field amplitude information, the experimenter 
can still retrieve the flash location and time of occur- 
rence; that is, (12) becomes a (3 X 3) system of linear 
equations, and f = col{x, y, d { ). 

If 8 is unity, the row vectors of K involving sin <j> t 
and cos</>, appear numerically small, that is, like a zero 


vector, relative to the other row vectors of K, and the 
matrix is ill-conditioned for many source locations when 
only three ALDF sensors are available. 

To avoid unstable inversions associated with an ill- 
conditioned K-matrix, we have made the assignment 8 
= 10 3 m. This increases the magnitude of the small 
trigonometric components of K and effectively filters 
small eigenvalues; see section 3d below and appendix 
A for additional details regarding the value of 8. Other, 
more sophisticated means of filtering small eigenvalues 
by adding external physical constraints to the solution 
process are discussed in Twomey (1977, Chapter 6). 


3. Simulated tests of the LP method 

a. Overview 

Because the effects of propagation are known to de- 
grade the quality of field amplitude data, F t [Cooray 
1987, his Eq. (2)], and because we have not taken spe- 
cific measures to model and correct for these errors, we 
will not consider these data in this and all tests to follow. 
We also assume that all sources and sensors are located 
on the surface of a spherical earth. By selecting a known 
source lat-long location, we generate the true arrival 
times and bearings to each sensor. Simulated measure- 
ments are generated by adding errors to the computed 
arrival times and bearings. The errors are chosen from 
a uniform random distribution. Although the bulk mag- 
nitude of errors that we will choose are typical of a real 
experiment, we have not attempted to simulate the de- 
tailed effects of propagation error. Such effects (with 
regard to arrival time retrievals) are provided by Honma 
et al. (1998). 

Next, the simulated measurements are analyzed with 
the LP method. Since the LP method is a planar model, 
we must establish a convention for mapping source and 
sensor locations (expressed in degrees of latitude and 
longitude) to locations in the x-y plane of a standard 
Cartesian coordinate system. We then apply the LP 
method to solve the problem in the Cartesian system. 
Next, an inverse mapping is used to convert the (x, y) 
solution back into latitude and longitude coordinates on 
the surface of the earth. At this point, the lat-long so- 
lution can be compared with the known source to assess 
true location error. 

If one assumes a flat earth and performs the entire 
simulation within a Cartesian coordinate system, the 
resulting retrieval errors are smaller. This is because one 
avoids errors due to earth curvature and the numerical 
errors associated with spherical/Cartesian system map- 
pings. Because in any real field experiment the source 
retrievals are ultimately referenced to the spherical 
earth, we include the net effects of earth sphericity in 
this and all other simulations provided below. 
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b. Spherical arrival time and bearing 


Figure 2 indicates how to compute the arrival time t t 
and bearing (/>, for the /th sensor on a sphere. The unit 
vectors pointing from the origin O to the /th sensor (M), 
to the lightning source (L), and to the North Pole ( N ) 
are, respectively, 

r, — cos<p, cosA,u + cos <p t sinA,v 4- sin(p ( w, 
r = cos<p cosAu 4- cos<p sinAv 4- sin<pw, 
f,v = w. (14) 

Using the law of cosines from spherical trigonometry 
gives the spherical angle A,-, 


A; = 


180 [cosa — cos b. cosc, 

COS 1 i 

it ^ sin b, sine, 


(15) 


where a = cos _1 (f N *f), b t = cos ! (f jV -f,), c t = 
cos _1 (r ( • r) = (1 IR)cti, R = radius of earth, and by 
convention the lightning source activates at t = 0. The 
angle, A,, varies between 0° (north) and 180° (south). 
We correct A, to construct the bearing function </>, that 
varies in the following manner: 0° (east), 90° (north), 
180° (west), 270° (south); that is, 


<p t - 90 - A, 
cf>, = 90 4- A, 
<j>, = 450 - A ( 


(northeast sources), 

(northwest and southwest sources), 
(southeast sources). (16) 


c. Mappings 

In general, different mappings produce different re- 
trieval errors. We consider two possible approaches: 
Mapping #1 (chosen for its mathematical simplicity), 
and Mapping #2 (chosen for its orthogonality). In Map- 
ping #1, we have 

x = (\ - A ,)R cos<p, y = {<P ~ <P\)R (17) 

where (A, <p) is an arbitrary longitude and latitude, re- 
spectively. The origin of the Cartesian coordinate sys- 
tem has been arbitrarily selected as site / = 1; that is, 
the ordered pair (A,, <p } ) is the location of site 1 andx(A 
= Aj, <p = <p { ) = 0, y(A = A,, <p = (pt) = 0. Note that 
y is measured along a great circle, that is, a longitude 
belt, but that x is measured along a latitude belt [which 
is only a great circle if <p x = 0 (the equator)]. 

In the second approach, or Mapping #2, we insist that 
both x and y are measured along great circles. To do 
this, we consider an orthogonal system (u, t, w) where 
u is a unit vector directed from the center of the earth 
to the intersection of the prime meridian and equator, 
w is a unit vector directed from the center of the earth 
to the North Pole, and v completes the ordered triple in 
accordance with the right-hand rule, that is, v = w X 
u. We then rotate this coordinate system through two 
Euler angles (A,, <p,) and define the new resultant 
(“starred”) system as (u*, v*, w*). In the starred sys- 
tem, u* is directed from the center of the earth to site 
1 . Mapping #2 is then 

x = /U*(A, <p) y = R<p*( A, (pi (18) 

where 


A*(A, (p ) = tan 


cos <p sinA cosA, — cos^ cosA sinA, 


^cos < p cosA cosA t cos<p l 4- cos <p sinA sinA, cos <p ] 4- sirup sincp 
<p*(A, <p) — sin ! (sin^ cos<p, — cos (p cosA cosA, sin<p, — cos <p sinA sinA, sin<p,). 


( 19 ) 


Again, one can verify from (18) and (19) that x(\ = 
A,, <p — <pj = 0, y(A = A,, <p = <p,) = 0. The arctangent 
expression in (19) must be appropriately corrected de- 
pending on what quadrant (northeast, northwest, south- 
east, southwest) the point (A, <p) is relative to site 1. 

d. Simulation 

We first consider three ALDF sites in the Darwin, 
Australia, region that were used as part of the MCTEX 
described in Keenan et al. (1994, 1996). Computer-gen- 
erated lightning sources were spaced 0.02° (~2 km) 
apart across the analysis region. Figure 3a shows the 


spatial distribution of the retrieved horizontal source 
location error, contoured in units of kilometers, when 
no experimental errors are considered and when Map- 
ping #1 is used. The retrieved location errors are within 
1 km for regions inside the ALDF network. 

Since no experimental errors have been added to the 
simulated values of the arrival times and bearings, the 
retrieval errors shown in Fig. 3a are due solely to earth 
curvature and numerical truncation error. We originally 
performed these simulations, and all simulations to fol- 
low, assuming a flat earth. When this was done, all of 
our methods gave retrieval errors well below 2.5 m 
across the entire analysis region when no measurement 
errors were involved. This is to be expected since our 
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Fig. 2. Spherical trigonometry used for determining arrival time 
and bearing. 

methods are exact solutions for the plane (the 2.5-m 
error maximum occurred only in the LP method over a 
limited portion of the analysis region and was an artifact 
of what accuracy level we required of our iterative ma- 
trix inversion routine). Hence, compared to the negli- 
gible errors obtained from the flat earth simulations, the 
errors shown in Fig. 3a are effectively due to earth cur- 
vature alone. However, the amount of retrieval error due 
to earth curvature depends on what cartesian-to-spher- 
ical coordinate system mapping is used (e.g., Mapping 
#1 or Mapping #2). Because one will always be inter- 
ested in how much retrieval error the planar models 
acquire due to earth curvature, all simulations below 
show retrievals first without added measurement errors, 
as in Fig. 3a. 

When experimental errors are included in the simu- 
lation, we obtain the retrieved location errors given in 
Fig. 3b. The retrieved errors are mean values obtained 
from performing 100 individual retrievals at each trial 
location. For each of the 100 trials, an arrival time error 
selected from a uniform random distribution (ranging 
from ~300 to 300 ns) is added to the arrival time value, 
and a bearing error (ranging from —2° to 2°) is added 
to the bearing value. In addition, we have simulated 
sensor location errors by purposely entering into the LP 
method false site locations (with an error as great as Vi 
m); the sensor location errors have remained fixed for 
all source analyses. As expected, the addition of ex- 
perimental errors increases location retrieval errors, but 
the retrieved errors are still within 10 km for a large 
portion of the analysis region. Roughly speaking (that 
is, not accounting for earth curvature errors, truncation 
errors, or other errors due to matrix inversion), a 300- 
ns timing error multiplied by the speed of light gives 
only a 90-m error, and a 2° error at a range of 300 km 
is about 10 km. 

When Mapping #2 is used instead of Mapping #1, 



Fig. 3. Lightning location retrieval errors for a three-station net- 
work using the LP method: (a) no measurement errors, and (b) with 
the following random measurement errors: 0.5-m sensor location er- 
ror, 300-ns timing error, 2° bearing error. Mapping is used. Con- 
tours arc in units of kilometers. The analysis region shown is 6° lal 
(667 km) X 6° long (about 651 km), and this is where the Maritime 
Continent Thunderstorm Experiment (MCTEX) was conducted. 


we obtain the results shown in Figs. 4a, b. As before, 
no experimental errors have been added to the sensor 
positions, arrival times, and bearings in the results of 
Fig. 4a, but the results in Fig. 4b include these errors. 
The results in Fig. 4a appear somewhat better than those 
in Fig. 3a, but the results in Figs. 3b and 4b are similar 
since the simulated experimental errors tend to mask 
differences between Mapping #1 and Mapping #2. 

Figures 5a-b and 6a-b show all of the same type of 
analyses just described, but for the case of four ALDF 
sensors. The additional sensor clearly helps reduce re- 
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Fig. 4. Same as in Fig. 3 except that Mapping #2 is used. 


trieval error. In addition, Mapping #2 produces smaller 
retrieval errors than Mapping #1. 

In summary, our simulations show that there are two 
basic principles operating here: 1) Mapping #2 generally 
outperforms Mapping #1, and 2) four-station networks 
are more resilient to measurement errors than three- 
station networks. Principle 1 is supported by the fol- 
lowing: (i) the error in Fig. 4a is less than in Fig. 3a, 
and (ii) the error in Fig. 6a is less than in Fig. 5a. Support 
for principle 2 is also evident. For the three-station net- 
work there are large increases in error going from Fig. 
3a to Fig. 3b and from Fig. 4a to Fig. 4b, but for the 
four-station network there is not as much increase in 
error (see error changes going from Fig. 5a to Fig. 5b, 
and from Fig. 6a to Fig. 6b). 

Next, a word regarding bearing data. When four sen- 
sors are available, bearing data are not needed to obtain 



longitude (degrees) 

Fig. 5. Same as in Fig. 3 except that four sensors are used. 

lightning location retrievals [see (11)]. By removing the 
bearing data from the four-station MCTEX region sim- 
ulations, that is assigning 5 = 0, and by applying Map- 
ping #2 we obtained virtually the same results as those 
presented in Figs. 6a,b. This is because the weighting 
factor 5 introduced into (11) to generate Figs. 6a, b is 
relatively small (i.e., 10*) so that bearing data has little 
influence on the final solution. We also find little change 
in the solutions for intermediate values, 5 = 10, 10 2 . 
However, as we increase 5 from 10* to 10 4 , 10 5 , and 
10 6 , the retrieval errors increase (see appendix A for 
more details). 

For a three-sensor ALDF network, bearing data plays 
a more profound role in helping to constrain the solution 
space. If 5 = 0 (no bearing data used), there would be 
fewer constraint equations than unknowns, and one 
would not be able to obtain a solution using the LP 
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longitude (degrees) 



longitude (degrees) 

Fig. 6. Same as in Fig. 3 except that four sensors and Mapping #2 
arc used. 


formalism. (Note: a different formalism to be described 
in section 4 below can be used to find solutions over a 
substantial portion of the analysis region using just three 
arrival time sensors.) As noted above, when 8 = 1 the 
K matrix is ill-conditioned for many source locations, 
leading to poor retrieval results. When an iterative meth- 
od is used to invert an ill-conditioned K-matrix, the 
computer time will be excessive. If a standard analytic 
form for K 1 is used (this is practical since the dimen- 
sions of K are small), then the time required to obtain 
the inverse is comparatively short and does not change. 
However, in either case, error magnification is exces- 
sive. When 8 = 10, KE, 10\ or 1 0 4 , there is no ap- 
preciable change in retrieval error and error magnifi- 
cation is minimal. When 8 = 10 s , the retrieval errors 
begin to increase slightly because of bearing errors. 



longitude (degrees) 



143.2 146.2 149.2 152.2 155.2 158.2 161.2 


longitude (degrees) 


Key: Location Error 



0-1 km 

EH 

1-5 km 

CU 

5-10 km 

I — 1 

>10 km 


Fig. 7. Same as in Fig. 3 except that this is for the TOGA COARE 
analysis region, and Mapping #2 is used. The analysis region is 18° 
lat (2002 km) X 18° long (about 1996 km). Shading, rather than 
contouring, is used to clarify the nonmonotonic distribution of re- 
trieval errors. 


Three-station LP simulations for different values of 8 
are also provided in appendix A. 

For comparison, we also provide error results (Figs. 
7a,b) for the three sites used in the Tropical Ocean and 
Global Atmosphere Coupled Ocean Atmosphere Re- 
sponse Experiment (TOGA COARE) described in Pe- 
terson et al. (1996) and Orville et al. (1997). This ex- 
periment employed a larger sensor baseline than that 
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used in MCTEX, and our simulated tests cover an anal- 
ysis region 18° X 18° in latitude and longitude. The 
sources in this simulation where placed 0.05° apart and 
Mapping #2 was used; as with the three-sensor MCTEX 
study, the retrieval errors for Mapping #2 differ little 
from those errors obtained using Mapping #1. Overall, 
the 2° bearing error and the effects of earth curvature 
make it difficult to obtain errors below 10 km for distant 
sources. In addition, the TOGA COARE network ge- 
ometry is not optimum since the baseline between sites 
1 and 2 is not very large in comparison to the other two 
baseline distances. 

4. Quadratic Planar (QP) method 

In this section we assume that only three arrival time 
measurements are available from the ALDF network. 
Hence, the methods of section 2 cannot be applied, but 
some insight about the source location can still be ob- 
tained. We assume that sensor / — 1 is at the origin of 
a rectangular Cartesian coordinate system (i.e., x x — y, 
= Z] = 0), we specify the convention f, = 0, and we 
again take z = 0. This leads to one nonlinear equation 
and two linear equations all in the three unknowns 

(x, }'> r): 

x 2 + y 2 = r 2 \ i - 1 

q\ = ~(r } - c 2 t~) 

= XjX + y t y + ct ( r\ i = 2, 3. (20) 

The equations in (20) were derived from the transit 
equation in (1) by means similar to that discussed prior 
to (2) of section 2. We have removed the source acti- 
vation time t with the relation t = —r/c. Since q ss 0, 
it is consistent that t < 0 given that r > 0. Each sensor 
is a distance r, = (x~ + yf) in from the origin. 

With minimal algebra, it can be shown that each equa- 
tion in (20) can be transformed into the equation of a 
circle with radius c(t, — t) and center at (x,, y ( ), where 
i = 1, 2, 3. We will see below that the source is located 
where the three circles intersect. We will also find that 
certain source locations produce arrival time data that 
can be described geometrically by two possible sets of 
three circles. Each set of three circles define a unique 
intersection point in the x-y plane, thereby leading to a 
fundamental ambiguity in source location retrieval. 

Geometric intersections of the circular curves de- 
scribed above are obtained by solving the system of 
equations in (20). To solve the system, we first subtract 
the terms ct t r from each side of the linear equation set: 

<7, , - q'i - ct,r = x ( x + y,y; i = 2, 3. (21) 

Identifying the vectors, q = col)#,, <7,), r = col(x, y), 
we may write 

q = Qr, (22) 

where the Q-matrix and its inverse are given by 


X 2 


Q - 1 

y> -y } z 

*3 

y* t 

“ > ? 2 

-x 3 x 2 


From (22) and (23) and our discussion preceding (21) 
we have the relations 

x(r) = O'., <?, - >2<?.0/O 2 y, - y 2 .v 3 ) 

y(r) = (x 2 q 2 - x i q t )/(x 2 y s - J 2 .v 3 ) 

t(r) = -r/c. (24) 

The x and y variables are written as functions of r in 
(24) since the components of q depend on r as given 
in (21). Substituting the first two equations of (24) into 
the first (nonlinear) equation of (20) and carrying out 
the algebra leads to an equation quadratic in r alone: 

Ar 2 + Br + C = 0, (25) 

where 


A = c 2 [r]t\ - 2(x 2 x 3 + y 2 y,)t 2 t, + r\t]] - (x 2 y, - y 2 x 3 ) 2 

B = 2c[—r\q\t 1 + (x 2 x 3 + y 1 y } ){q\t 3 + q{t 2 ) ~ rlq\t 3 ] 

C = r 2 q[ 2 - 2(x 2 x 3 + y 2 y,)q r 2 q\ + r 2 q’ 2 . (26) 

Hence, the lightning source range, r, is the nonnegative 
real root obtained from the formal (two root) solution: 




r = 


r 


-B + VS 2 - 4 AC 


2A 

-B - VB 2 - 4AC 
2 A 


(27) 


Values of r = 0 correspond to a direct lightning strike 
of sensor i = 1, which we ignore. Note from (26) that 
the numerical value of the coefficients (A, B, C ) are 
obtained from the sensor locations and excitation times, 
that is, on the six variables: {x 2 , y 2 , t 2 , x 3 , y u r 3 ); the 
variables, q\, are obtained from the expressions l(r? — 
c 2 tj) as given in the last two equations of (20). After 
these data are used to compute r, (24) is used to find 
the lightning location [x(r), y(r)], and time of occur- 
rence, t(r ). 


5. Simulated tests of the QP method 

By placing computer-generated lightning sources 
0.02° apart in latitude and longitude across the analysis 
area (see section 3), we have determined the horizontal 
location error resulting from each root in (27). To fa- 
cilitate comparisons with simulated tests of the LP meth- 
od, sensor position and arrival time errors used here are 
as described in section 3d and 100 trials at each source 
location are once again used to generate mean retrieval 
location errors. Mapping #1 given in (17) is employed. 

Figure 8 clarifies what root provides the correct 
source location. The shaded regions are where r + is 
correct and the unshaded region is where is correct. 
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Interestingly, the dividing lines of these regions are de- 
fined by the sensor baselines. 

When we pick only the correct root and plot the as- 
sociated error result over the analysis region, we obtain 
the result given in Fig. 9a. When a 300-ns uniform 
random error is added to the computer-generated arrival 
times, we obtain the mean horizontal distance errors ~ 

given in Fig. 9b. Considering that only three sensors E 

are involved, retrieval errors are quite good; a large 3 

region of errors below 1 km is evident. Distant sources, J 

or sources located near the outer sensor baselines are 
more difficult to accurately retrieve. See appendix B for 
a definition of outer sensor baseline . A comparison be- 
tween the three-station LP results in Fig. 3b shows that 
the QP method provides better results over most of the 
analysis region. This is due in part to the large 2° bearing 
errors implemented in the LP simulation and the fact 
that the three-station LP method depends on bearing 
data to obtain a solution. The four-station LP method 
(Fig. 5b) outperforms the QP method in most locations. 

Since bearing data can aid in determining which root, 
r + or r_, is correct (see section 6a below on solution 
ambiguity) and since the QP method generally gives 
better results than the three-station LP method, it is 
evidently better to use the QP method than the three- 
station LP method when bearing data are available (pro- 
vided the source is not located near an outer sensor 
baseline). This conclusion is based, of course, on an 
assumed bearing error of 2°. 

6. Examination of QP method roots 

When applying the QP method to actual arrival time 
data, one picks the solution associated with a nonneg- 
ative real root, that is, the source range r must be non- 
negative and real. A detailed discussion of root results 
is provided below. 

a. Unequal nonnegative real roots (ambiguities) 

In the discussion of three sensor networks by Holle 
and Lopez (1993, pp. 8, 11), the lightning source lo- 
cation (x, y) is described in terms of the intersection of 
two hyperbola branches; each branch is defined by two 
sensors. For some lightning source locations, the hy- 
perbola branches intersect at two locations [see, for in- 
stance, Fig. 6, p. 11 of Holle and Lopez (1993)]. This 
amounts to a fundamental ambiguity in location retrieval 
and the authors correctly assert that the ambiguity can 
be removed by adding a fourth (properly positioned) 
sensor. Ambiguities are described in our formalism by 
the intersection of circles as indicated in section 4. To 
fully appreciate the hyperbolic and circular geometrical 
viewpoints, it is important to recognize that the two 
intersection points defined by two sets of three circles 
are identical to the double intersections obtained from 
the two hyperbola branches mentioned above. In other 



I 1 : r + correct 

I 1 : r_ correct 


Fig. 8. Comparison of roots in the QP method. The shaded regions 
indicate where root r _ is the correct solution and the unshaded region 
indicate where root r is the correct solution. No measurement errors 
have been added to the simulated data. 

words, these two widely different geometrical view- 
points produce identical results, as they must. 

Without reference to the geometry of hyperbolic or 
circular intersections, our algebraic formalism imme- 
diately defines all ambiguous cases. An ambiguity will 
exist whenever two unequal nonnegative real roots re- 
sult from (27). In order to determine what lightning 
source locations produce these “ambiguity regions,” we 
have kept a record of the root results in the numerical 
experiments described in section 5 above. For the case 
of no simulated experimental errors, the source locations 
that resulted in two unequal nonnegative real roots are 
indicated by the shaded regions in Fig. 10; see section 
6b below for minor corrections to the ambiguity regions. 
In general, a different network geometry would produce 
different results. 

Strictly speaking, since two distinct sources can pro- 
duce identical arrival time difference information, there 
is no means of discriminating which source location is 
correct unless some additional information is supplied 
to the retrieval process. In effect, the solution is fun- 
damentally nonunique. [Similar comments about the re- 
trieval of charge from ground-based field measurements 
have been made in Koshak and Krider (1994). In that 
problem, a point charge Q u and a sphere of radius a 
with total charge Q t> produce identical electrostatic fields 
outside the radius a.] Hence, additional measurements 
(e.g., arrival time, bearing, signal amplitude, radar, 
acoustical, interferometric) must be used to pick the 
correct root. Bearing data would be the most common 
data to use in root discrimination since it is part of the 
ALDF data stream. 
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longitude (degrees) longitude (degrees) 



longitude (degrees) 

Fig. 9. Retrieval errors from the correct root in the QP method: 
(a) no measurement errors, and (b) with the following random mea- 
surement errors: 0.5-m sensor location error, 300-ns timing error. 
Mapping #1 is used. Contours in kilometers. 


Nonetheless, an experimenter might be tempted to 
compare the shaded regions in Figs, 8 and 10 in order 
to determine which of the two unequal nonnegative real 
roots produce the true source location. However, one 
must remember that Fig. 8 does not pose a real physical 
constraint to an unknown source since it only provides 
information if the source location is already known ; 
obviously in a real field experiment the source location 
is not yet known. 

Additional rigor clarifies the immutability of the am- 
biguous case. Note that there are three ambiguity regions 
in Fig. 10, and call the regions {R ,, R z , P 3 }. Similarly, 
there are three regions {P ] , P 2 , P,} in Fig. 8 that are 
wholly contained within each of the respective ambi- 


I I : ambiguous 
I I : not ambiguous 

Fig. 10. Shaded regions are the QP method “ambiguity regions” 
that indicate what lightning source locations result in two unequal 
nonnegative real roots. No errors were added to the simulated arrival 
times. 

guity regions. If one subtracts the respective regions 
{P,, P 2 , P 3 ) from the respective regions {P M P 2 , P.}, 
one obtains the three regions {N ly N 2 , N z } = {P, — 
P, , R 2 - P 2 , P 3 — P 3 } . When two unequal nonnegative 
real roots are obtained, we find that r + occurs in P t and 
r_ occurs in N t , that is, the pair of roots produce so- 
lutions in (P,, JV,), (P 2 , AP), or (P 3 , AT 3 ). From Fig. 8, 
region P. is the region where r+ is correct, N t is the 
region where r_ is correct, and both P. and N t are sub- 
regions of R r In other words, each solution is possibly 
correct so that comparisons between Figs. 8 and 10 serve 
no help in determining the correct root. 

Nonetheless, r, is correct for most of the ambiguity 
region shown in Fig. 10. Hence, a MCTEX experimenter 
who obtains two nonnegative real roots, but does not 
have ancillary datasets such as radar, magnetic bearing, 
etc. (for determining the correct root) is best off se- 
lecting the root r 

b. Equal nonnegative real roots 

In this section we are interested in identifying what 
source locations produce two equal non negative real 
roots. Note that this condition is satisfied when the dis- 
criminant, B 2 — 4 AC, in (27) is zero, that is, the two 
equivalent roots correspond to a unique (unambiguous) 
solution [.v(r), y(r), t(r)] where r t = r_ — r. In appendix 
B, we show for arbitrary noncollinear network geom- 
etries that the discriminant function is zero only along 
the outer sensor baselines. Therefore, the ambiguity re- 
gions shown in Fig. 10 are technically not ambiguous 
along these linear domains. 
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Fig. 11. Plot of the QP method discriminant (scaled by a large 
arbitrary constant for plotting purposes). Expressions for A, B , and 
C in the discriminant, B 2 - 4 AC, are given in (26). 

c. Complex roots 

Complex roots occur whenever the discriminant in 
(27) becomes negative. Figure 1 1 shows how the dis- 
criminant varies for different source locations across the 
analysis region. As we have already shown in section 
6b and appendix B, the discriminant is zero for sources 
located along the outer sensor baselines. Figure 1 1 
shows additionally that the discriminant is a relative 
minimum at the outer sensor baselines. 

From the simulation in section 5 (with 300-ns arrival 
time errors, and 100 trials per test location) we have 
tallied the fraction of trials at each location that produce 
complex roots. Figure 12 shows that there are no com- 
plex roots over most of the analysis region except when 
the sources are near the outer sensor baselines. These 
regions (or “spokes”) appear to diverge with range from 
the sensors and as many as 40%-6O% of the sources 
are complex within the spokes. Clearly, for sources lo- 
cated sufficiently close to the outer sensor baselines, 
measurement errors are occasionally large enough to 
drive the discriminant negative. Whenever the discrim- 
inant is negative, both roots in (27) are complex and no 
physical solution is obtained. Conversely, whenever 
complex roots are obtained from a set of actual mea- 
surements the source is likely to be located in one of 
the spoked regions. 

d. Overview of root results 

From our discussion so far, we can conclude that any 
retrieval will produce one of the following cases: (a) r r 
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Fig. 12. Fraction of 100 simulated sources at each location that 
produce complex roots using the QP method. A 0.5-m sensor location 
error, a 300-ns timing error, and Mapping #1 was used. 


> 0, r ^ 0, r, 7 ^ r_; (b) r i > 0, r_ > 0, = r ; 

(c) r + < 0, r_ ^ 0; or (d) r + and r_ complex. Case a 
corresponds to a source that is located inside one of the 
ambiguity regions, case b corresponds to a source lo- 
cated on an outer sensor baseline, case c corresponds 
to a source that is not located in any of the ambiguity 
regions or along any outer sensor baseline, and case d 
corresponds to a source located on or near any outer 
sensor baseline when measurement errors are sufficient 
to drive the discriminant negative. 

Note that we do not include the case r + > 0, r_ < 
0 since if r + > 0, the source must lie in one of the 
ambiguity regions implying that would be nonneg- 
ative (i.e., a contradiction). We also disregard the case 
that both roots are negative since a physical source must 
lie a nonnegative distance from sensor / = 1. 

7. Sample storm analyses 

We have applied the LP and QP algorithms (with 
Mapping # 2 ) to retrieve several thousand cloud-to- 
ground flashes that occurred over the MCTEX analysis 
region during 28 November (day 332) and 29 November 
(day 333) 1995. To demonstrate the internal consistency 
between various forms of the algorithms and to clarify 
to what degree the algorithms differ, we have analyzed 
only those flashes that produced a valid excitation at all 
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Fig. 13. Day 332 ground flash retrievals using the following algorithms: (a) LP2, (b) LP3, (c) 
QP, and (d) LP4. Only those source retrievals falling inside the above MCTEX region arc shown. 
Of the four algorithms, the results of LP4 are considered the best estimate of the true lightning 
locations. 


four of the ALDF sites. In this way, each flash can be 
analyzed by each selected algorithm. 

The LP two-station algorithm, or “LP2” algorithm 
ingests the arrival time and bearing data from sites 1 
and 2, and ignores the data from sites 3 and 4. The LP3 
algorithm ingests the arrival time and bearing data from 
sites I, 2, and 3, and ignores site 4. The QP algorithm 
ingests only arrival time data from sites 1, 2, and 3, and 
ignores site 4. Finally, the LP4 algorithm ingests arrival 
time and bearing data from all four sensors. Hence, the 
results derived from LP4 are considered to be the best 
retrieval results against which the three remaining al- 
gorithms are compared. For example, the discrepancy 
between LP2 and LP4 results indicates how appropriate 
it might be (in future storm analyses) to use LP2 when 
only two sites detect the flash. 

Figure 13 shows the ground flashes derived from day 
332 for each algorithm, and Fig. 14 shows the results 
derived from day 333. Note from the LP2 results pro- 
vided in Figs. 13a and 14a that there are obviously 
location errors along the baseline between site 1 and 2, 
in agreement with the general error results provided in 
section 3. The cluster of flashes to the south on day 332 
(Fig. 13d) is reasonably well retrieved by LP2 (Fig. 13a) 
since the cluster is not near the baseline, but it is some- 


what smeared out and biased closer to the network. Sim- 
ilarly, the myriad of flash clusters to the south on day 
333 (Fig. 14d) are reasonably well retrieved by LP2 
(Fig. 14a), but again there is positional smearing. 

The results of LP3 (Figs. 13b and 14b) dramatically 
reduce the baseline errors and smearing that was as- 
sociated with the LP2 results. However, there are still 
some positional adjustments between the LP3 and LP4 
results. 

Finally, the results of the QP algorithm are shown in 
Figs. 13c and 14c. In all of the QP results shown here, 
root ambiguities were resolved using the LP3 source 
solutions. Generally, these results agree favorably with 
the LP4 results and are perhaps even better correlated 
to the LP4 results than are the LP3 results, for some 
source locations. However, several sources incorrectly 
fall along an arc-like pattern that closely resembles the 
boundaries of the ambiguity region provided in Fig. 10. 
In the ambiguous case, the simulations of section 6a 
also showed us that one solution occurs near the bound- 
ary of the ambiguity region and the other well inside 
the (same) ambiguity region. Because of measurement/ 
baseline errors and our planar assumptions, LP3 did not 
always pick the correct root. Hence, sources occurring 
within the southern cluster (the solution associated with 
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Fig. 14. Day 333 ground flash retrievals using the following algorithms (see algorithm de- 
scriptions in text): (a) LP2, (b) LP3, (c) QP. and (d) LP4. Only those source retrievals falling 
inside the above MCTEX region are shown. Of the four algorithms, the results of LP4 arc 
considered the best estimate of the true lightning locations. 


the correct root) had retrieved locations near the bound- 
ary of the ambiguity region (the solution associated with 
the incorrect root). (Note in fact that the southern cluster 
of Fig. 13c consists of far fewer flashes than in Figs. 
1 3a,h,d indicating that there has been a shift of solutions 
from the cluster to the border.) In addition, the increase 
in retrieval errors near the outer sensor baselines (see 
Fig. 9b) are large and also play a role in defining how 
well the cluster position is retrieved. Finally, sources 
very near the outer sensor baselines created complex 
roots (no solutions). Because the LP4 and QP methods 
are significantly different mathematical approaches, it 
is particularly encouraging to see such good agreement 
between them. 

8. Summary 

In this writing we have derived, tested, and applied 
two basic methods for retrieving the location and time 
of occurrence of lightning ground strikes from a network 
of ALDF sensors. We developed the methods so that 
they could be used in future validation studies of NASA 
space-borne lightning sensors: the OTD and the LIS. 
Because these imagers have a nadir resolution of 8 and 
4 km, respectively, our retreival algorithms are adequate 


for validating lightning source locations within a few 
hundred kilometers of the ALDF network. As such, we 
have avoided the need to use the (proprietary) nonlinear 
X 2 minimization algorithm mentioned in Cummins et 
al. (1993, 1995, 1998). The development of our algo- 
rithms has also expanded upon and clarified the theo- 
retical aspects of ALDF data inversions. 

The first approach introduced in this writing, or LP 
method, assumes that arrival time, bearing, and field 
amplitude measurements are all available from the net- 
work. As provided in (12), these measurements are col- 
lected into one coherent linear system of equations that 
is solved by straightforward inversion. Because of the 
general form of (12), we have clarified what solution 
options one has if only a subset of the measurements is 
available (as when sensors do not trigger on a distant 
or low amplitude event and/or when sensors are defec- 
tive). When only three arrival time measurements are 
available, the LP method cannot be used. To solve this 
problem, we have introduced a second approach, or QP 
method. 

Both the LP and QP methods described in this writing 
have problems for sources located near or on the outer 
sensor baselines. In the case of bearing data, the bearing 
lines-of-site intersect at small angles, and this fact gives 
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rise to large location retrieval error. In the case of timing 
data, the algebraic solution space is defined geometri- 
cally in terms of the intersection of circles (or equiva- 
lently, as hyperbolic branches). For sources located near 
or on the outer sensor baseline, the circles (or hyperbolic 
branches) intersect at small angles, and again give rise 
to large location retrieval error. For three- and four- 
station networks, the most accurate solutions are found 
in the region bounded by the three outermost sites. For 
two-station networks, the most accurate solutions are 
found outside the line segment joining the two sensors. 
For all networks, location retrieval error increases with 
source range. 

Since a lightning flash need not trigger all sensors in 
a given network of n sensors, all of the solution tech- 
niques we have described need to be considered. For 
our four-sensor MCTEX network, we have found that 
it is best to use LP4 when all four sensors trigger on a 
discharge. If only three sensors trigger on a flash, we 
use the QP method (unless the flash is located near an 
outer sensor baseline, in which case we would employ 
the LP3 method). LP3 solutions would sometimes be 
favored over those QP solutions that align along char- 
acteristic “arc patterns” such as those shown in Figs. 
13c and 14c. Finally, we employ LP2 when only two 
sensors detect a flash; because this method has difficulty 
when the flash is located anywhere along a line passing 
through the two sensors, we will consider constraining 
the solution with signal amplitude and/or other ancillary 
datasets. 

The planar methods express the source locations di- 
rectly in terms of the measurements. The solutions are 
concise, require little computer time, and afford the user 
with specific physical insights about the retrieval prob- 
lem in the form of analytic equations. The equations 
express the relative importance and effects of timing/ 
bearing data on the accuracy of final solutions. They 
also describe the regions of ambiguity, and the regions 
where source locations produce complex roots. 

Moreover, the LP and QP methods introduced here 
offer the authors and other researchers a means to in- 
tensively analyze and compare, first hand, lightning- 
radio-source locations with OTD-LIS low-earth orbit 
lightning detections. In the future, we intend to apply 
these methods to analyze a wide range of thunderstorms, 
to continue intercomparing the methods, and to relate 
the results to OTD and LIS and other independent da- 
tasets such as radar, Lightning Detection and Ranging 
(LDAR), and the National Lightning Detection Network 
(NLDN). The first author will also improve some of the 
matrix methods presented here to directly account for 
earth sphericity; more elegant oblate spheroidal models 
are also under consideration. 
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APPENDIX A 

Weighting of Bearing Data 

It is of interest to determine to what extent bearing 
data is actually being used to constrain the lightning 
source location in the LP method as a function of the 
weighting factor 8 introduced into the linear system of 
equations provided in (11). Insight is gained by con- 
sidering the case of three sensors and ignoring field 
measurements, F t . Then (II) reduces to 


«2 


(* 2 ~*i) ( y 2 — i ) c(r, - / 2 ) 
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This system, which can be written in the standard notation 
g = Kf, has the least squares solution f = (KK) 1 Kg. We 
are interested in the explicit functional dependence of x 
and y on the arrival time t, and bearing (ft , data, and on 
the weighting factor 8. Because this is a very involved 
hand calculation, we utilize a computer-aided symbolic 
manipulator to arrive at the following form: 


x = 


1 

h 0 8 2 + h x 


h 2 a 2 + /7 3 a 3 + 2 (M 2 + 

j-4 


(A2) 


The 10 functions, hj (j = 0, . . . , 9) depend, in general, 
on the arrival time and bearing data. The variables (a 2 , 
a 3 ) depend only on arrival time data and network ge- 
ometry, and the variables (j3,, /3 2 , jS.O depend only on 
bearing data and network geometry as given in (3) and 
(6), respectively. A similar form holds for y. The co- 
efficient in front of the square brackets in (A2) does not 
preferentially weight the a’s or /3’s so it is of no concern 
in this discussion. However, the coefficients (/? y + 
8 2 h J+i ) do weight the /3's but not the a’s (i.e., these 
coefficients weight the bearing data, but not the arrival 
time data). This leads to the final results 
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longitude (degrees) 


8 = 0 => no solution 
lim .v = -j-( hiP, + h s P 2 + = A($,); 

Rq 

i = I, 2, 3. (A3) 

The first result is true since when 6 = 0 there are two 
equations in three unknowns in (Al). The second result 
is true because the only arrival time dependence asso- 
ciated with each function ( h 0 , /? 4 , h 5 > h b ) is a factor 
{t\ + *3)- When the ratios hJh Q are taken, this factor 
cancels out and we are left with a function, A, that 
depends only on bearing data. 

In summary, for a positive finite value of 8, both 
arrival time and bearing data are utilized. How r ever, as 
8 is increased from zero, bearing data eventually be- 
comes more heavily weighted over arrival time data 
until, for a sufficiently large value of 6, only bearing 
data is being used to determine the source location. For 
8 = 1, we have difficulty inverting K for many source 
locations. We have also performed retrievals, in the pres- 
ence of measurement errors, for the values: 8 = 10, 10 2 , 
10\ 1 0 4 , I0\ and 10 6 . Figure Al shows the result for 
8 = 10, 10\ and 10 5 . There is not much change in the 
solution from 10 to 10\ but the dominance of bearing 
data constraints at 10 5 and higher begins to reduce the 
quality of the solution (i.e., a 2° random bearing error 
can create a substantial location error if the source range 
is sufficiently large). 

We have performed the same type of computer-aided 
symbolic manipulation to determine explicit forms 
when a four-sensor network is used [i.e., one more ar- 
rival time equation and one more bearing equation are 
added to the system in (Al) so that K becomes a 7 X 
3 matrix]. In this case, the form of * (as well as y) is 


1 

k 0 8 4 + M 2 + K 


x 


2 (kj+i 8 2 + kj+4) a j + 2 + kj+\ 2 8 2 )Pj 

;=2 /=i 


(A4) 


Most of the 17 functions, k J9 j = 0, . . . , 16 depend on 
both arrival time and bearing data. However, (k 2f k b , k 7 , 
k^) depend on arrival time data, but do not depend on 
bearing data. We obtain the final limiting conditions 


i— 


Fig. AL Location retrieval errors for a three-station network using 
the LP method when (a) 8 = 10, (b) 5 = I0\ and (c) 5 = I0\ The 
same simulated measurement errors discussed in section 3 are used: 
0.5-m sensor location error, 300-ns timing error, 2° bearing error. 
Mapping #2 is used. Contours (km). 
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lim x - -j~(k b a 2 + k 7 a 3 + k H a 4 ) = r(f ( -) 

5— >0 ki 

lim .v = — (fc.,/3, + k m /3 2 + fc n j3, 4- ^ l2 jSJ = ft (<£,); 

5^-jr Ky 

r = 1,2, 3. (A5) 

Hence, we swing from a solution governed only by 
timing data (5 = 0) to one governed only by bearing 
data (8 = & ) . The solutions for 8 = 10, 10\ and 10 5 
are shown in Fig. A2. The bearing data significantly 
worsens the solution when weighted heavily (8 = 10 5 ). 

APPENDIX B 

Locations Where QP Method Discriminant 
Function Vanishes 

We investigate more rigorously the zeroes of the dis- 
criminant function A = (B 2 — 4 AC) of (27) in the QP 
method. The computer plots of this function gave some 
interesting results near the outer baselines of the sensors, 
that is, there appears to be minima there (see Fig. 11). 

In the following formal approach, we algebraically 
reduce the discriminant into the product of three factors. 
Each factor is then shown to vanish along a specific 
outer sensor baseline. Our results apply to arbitrary net- 
work geometries. Because a zero discriminant implies 
that two nonnegative, real, and equal roots are obtained, 
a unique (unambiguous) solution [*(r), y(r), r(r)] is ob- 
tained on the outer sensor baselines where r = ~BI(2A) 
= r+ = r_. 

Using the forms in (26) for A, B, and C, the discrim- 
inant can be written as 

A = 4c'-[(p'- - rlr^lq'rt] - 2 q’ : q[t 2 t } + q?tj) 

+ e 2 &], (Bl) 

where 

P = r 2 • + 3^3 

1 _ 1 

s = - detQ = -(* 2 y 3 - y>* 3 ) 
c c 

<7 = rlq[ 2 - 2pq\q\ + r\q \ 2 , (B2) 

To simplify some of the algebra without losing gener- 
ality, we rotate the x and y axes so that y 2 = 0. Further 
reduction of (Bl) leads to 

A = [*Iy, J C*i - c 2 tl)][r 2 3 ~ c 2 f|] 

X [(xl - c 2 t\ ) + 2(c 2 /,fj - + (r| - 

(B3) 

Figure Bl considers the second factor, (rf - c 2 /J), in 
(B3). For a source located on the solid line with r s r 3 
site three is excited at 




longitude (degrees) 

Fig. A2. Same as in Fig. Al, but for a four-station network. 
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A 

y 



Fig. B I . Geometry for analyzing the second factor in the discrimi- 
nant function. 


A 


y 



Fig. B2. Geometry for analyzing the third factor in the discrimi- 
nant function. 



The factor becomes 


(B4) 

c 



Proceeding in a similar fashion, the factor for a source 
on the dashed line with 0 < r < r 3 is the nonzero result, 
4 r(r 3 — r), and the factor for a source on the dotted line 
(with 0 < r < r 3 ) or a source on the thick line (with r 
> rf) is zero. Hence, the factor is zero along the line 
running through the sensors (including the sensor lo- 
cations themselves, but excluding the line segment be- 
tween the sensors). This is what we refer to as the “outer 
sensor baselines. ” Similar comments can be made re- 
garding sites 1 and 2 when the factor (x| — c 2 t\) is 
considered. 

Evidently the third factor in (B3) corresponds to the 
line running through sites 2 and 3. To prove this, we 
consider the geometry provided in Figure B2. For a 
source on the solid line with d > 0 we have 


D + d — r 



(B6) 


Substituting these expressions into the third factor, and 
noting that D 2 = (.v 2 — jc 3 ) 2 + yj = x\ — 2x 2 x 3 + rj, 
we obtain 


x\ - (D + d - r) 2 + rf - (d - r) 2 
-F 2 (D + d — r)(d — r) - 2x 2 x 3 


= (xj - 2x 2 x 3 + r[) - [ (D + d - r) - (d - r)] 2 
= (xj - 2x 2 x 3 + rf) - D 2 = 0. (B7) 

For a source located on the dashed line (but not at site 


2 or site 3) the factor reduces to the nonzero result: 
4(x) 2 <x) } [(x 2 — x 3 ) 2 4- y 2 ], where the constant factors (w 2 , 
co 3 ) obey the constraints to 2 + co 3 = 1, u) 2 > 0, o> 3 > 
0. Finally, for a source on the dotted line and a distance 
I > 0 from site 2, we have 


h 



D + l - r 


(B8) 


These expressions have the same form as those in (B6) 
but are interchanged. When substituted into the third 
factor of (B3), the factor reduces to zero as in (B7). 
This completes the proof showing that the discriminant 
function vanishes along the outer sensor baselines. 
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